﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;

namespace ProjectEulerSolutions
{
    /*
     * Consider quadratic Diophantine equations of the form:

x^2 – Dy^2 = 1

For example, when D=13, the minimal solution in x is 649^2 – 13×180^2 = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

3^2 – 2×2^2 = 1
2^2 – 3×1^2 = 1
9^2 – 5×4^2 = 1
5^2 – 6×2^2 = 1
8^2 – 7×3^2 = 1

Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

     * */
    class Problem66 : IProblem
    {
        public string Calculate()
        {

            //wikipedija kaže da se to zove Pell's equation
            //rješenje je naći prvog konvergenta hi/ki od continued fration od sqrt(d) koji zadovoljava jednadzbu
            //hi^2 - d*ki^2 = 1

            //a ne ic zivom silom... nebi nikad, a ovako 70ms

            int maxD = 0;
            BigInteger maxX = 0;

            for (int d = 2; d <= 1000; d++)
            {
                BigInteger x = 1;
                BigInteger y = 1;

                if (Math.Sqrt(d) == (int)Math.Sqrt(d))
                    continue;

                int i = 0;

                while (x * x - d * y * y != 1)
                {
                    i++;
                    GetNthConvergent(d, i, out x, out y);
                }

                if (x > maxX)
                {
                    maxX = x;
                    maxD = d;
                    Console.WriteLine("x:{0}, d:{1}, y:{2}", x, d, y);
                }
            }

            return maxD.ToString();
        }

        void GetNthConvergent(int d, int n, out BigInteger h, out BigInteger k)
        {
            int index = n;

            //oooooogromni brojevi...
            BigInteger nom = 0;
            BigInteger denom = 1;

            List<int> cf = GetContinuedFractionForRoot(d).ToList();

            while (index >= 1)
            {
                int whole = 0;
                if (index > cf.Count)
                {
                    int target = (index - 1) % (cf.Count - 1);
                    whole = cf[target];
                }
                else
                {
                    whole = cf[index - 1];
                }

                BigInteger tempNom = nom;
                nom = denom;
                denom = whole * denom + tempNom;
                index--;
            }

            h = denom;
            k = nom;
        }

        IEnumerable<int> GetContinuedFractionForRoot(int n)
        {
            int nom = 0;
            int denom = 1;
            List<int> CFs = new List<int>(); //Continued Fractions
            double root = Math.Sqrt(n);

            if (root != (int)root)
            {
                int first = (int)((root + nom) / (double)denom);
                while (true)
                {
                    int whole = (int)((root + nom) / (double)denom);

                    //Console.WriteLine("{0} / {1}", nom, denom);

                    yield return whole;

                    int result = nom - whole * denom;
                    nom = -result;
                    denom = (n - result * result) / denom;

                    if (denom == 1)
                    {
                        yield return first + nom;
                        break;
                    }

                }
            }
        }
    }
}
